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ABSTRACT 


A new approach is proposed for the design of approximate, fixed-order, 
discrete time realizations of stochastic processes from the output covariance 
■)SCi,j) over a finite time interval. No restrictive assumptions are imposed 
on the process; it can be nonstationary and lead to a high dimension realiza- 
tion. Classes of fixed-order models are defined, having the joint covariance 
matrix of the combined vector of the outputs in the interval of definition 
greater or equal than the process covariance (the difference matrix is non- 
negative definite). The design is achieved by minimizing, in one of those 
classes, a measure of the approximation between the model and the process 
evaluated by the trace of the difference of the respective covariance matrices. 
The models belonging to these classes have the notable property that, under 
the same measurement system and estimator structure, the output estimation 
error covariance matrix computed on the model is an upper bound of the corre- 
sponding covariance on the real process. 

An application of the approach is illustrated by the modeling of random 
meteorelogical wind profiles from the statistical analysis of historical data. 


I. INTRODUCTION 


Given an m- vector- valued discrete-time process y(*) over a finite 
interval [1,N], with a positive definite eoVariance ^(i,j), we consider the 
problem of computing a markovian representation, with output y^Ci) . which 
statistically approximates y(i) in the interval 1 S i 5 N of the form 


x(i + 1) = A(i) • x(i) + B(i) • w(i) 
ynCi) = C(i) • x(i) + D(i) • w(i) 


( 1 . 1 ) 


where the dimensions of x(i), w(i), A(i), B(i), C(i), D(i)^ are nxl, pxl, 
nxm, nxp, mxm, mxp, respectively, and w(i) is a p-vector-valued independent 


^The matrices I>(i) can also be assumed to equal zero. 


noise. The dimension p of w(*) and the order n of the model are chosen 
by the designer, noting that the order of the model may be lower than the 
exact minimal realization of the process. Neither stationarity nor the 
hypothesis on the order of the minimal realization of the process is 
postulated. 

The model is used to design and evaluate linear recursive estimators of 
the process output in the interval 1 S i 5 N - that is, the process output 
will be observed through a linear measurement system with independent additive 
noise. The measurement system and, in turn, the filter characteristics are 
not actually specified. However, for the same arbitrary measurement system 
and estimator, the output estimation error covariance matrix^ computed on the 
model from measurements in the interval 1, . . . , i (1 5 i - ND is required 
to be greater than or equal to the equivalent estimation covariance on the 
real process, that is, the difference matrix is nonnegative. The models with 
this property are indicated here as "statistically guaranteed models" of the 
process. 

General solutions of the exact or partial realization (specifically 
minimal realization) problems exist for stationary processes with a lumped 
representation (refs. 1-5). In contrast, dealing with nonstationary processes, 
the literature commonly assumes that an analytic expression of the output 
covariance is available in a separable form (refs. 6-9). This form also 
implies a lumped realization. 

Therefore, if only the numerical values of the covariance are given, the 
computation of a separable form, when it exists, is preliminarily required. 

This approach, though theoretically viable, does not work out the practical 
difficulties encountered when the resulting realization is of high dimension 
or of changing structure in time. Dealing with an experimental output covari- 
ance which does not satisfy precise analytical assumptions, it may be necessary 
to adopt a lower-order approximation in several practical applications. How- 
ever it seems that no^ general results are available to the approximate 
realization of an arbitrary process with a constrained-order model. 

The approximation technique proposed here is to minimize in a class of 
"guaranteed models" of fixed order n a measure of the approximation between 
model and process evaluated by the trace of the difference between the respec- 
tive covariance matrices (a partial minimization scheme is proposed in the 
paper). It will be shown that guaranteed models are characterized by having 
an output covariance matrix greater than or equal to the process covariance.** 
This property offers, among other advantages, the option to use the relatively 
simple performance measure that we have proposed. It is evident that computa- 
tional limitations exist if N is large and no restrictive assumptions are 
taken for the process. It is shown, however, that the general approach can be 

^The joint covariance matrix of the combined vector of filtering and 
prediction output errors from i to N. 

^For an example of approximation of stationary processes by ARMA models, 
see reference 10. 

**A similar property in a less general approximation approach was already 
mentioned in reference 11. 
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applied with simplifications to stationary processes with finite dimension 
realization even if N is large. 

The probabilistic model of vertical profile of atmospheric wind velocity 
from statistics of historical records of wind measurements is a case of non- 
stationary, high-order realization^ process and it is presented here as an 
example of the method. The model has been used in refere^ice 12 to estimate 
the onboard wind shear estimation from an airplane in descent flight. 


II. GUARANTEED COVARIANCE ESTIMATION 


Consider the estimation at instant i (1 < i 5 N), with a linear 
estimator, of the present and future portion of the process y(j) (j = ij • • •> 
N) from the present and past measurements operated on the process output through 
a linear measurement system with additive independent noise. Referring to 
past portions of the process output, define the combined vector 


= [y'Ci), y'(i - 1), . 

and the joint covariance matrix 

"'«(i,i), . . . , KCi,l) 


. , y’CD] 


« 6 


C2.1) 


P(i) = 


.... 


1 < i < N 


(2.23 


with a similar notation defined for the model The present 

and future portion of the process is indicated as 

' ( 2.33 

where ^2 a matrix (N - i + 13m x Nm of the form 

(?2 = [^CN - i + IJmjO] (2.43 

and its estimate, by a linear estimator is 

= ^ • Z (2.53 

The vector Z of dimension m * i is the present anc past noisy measurement of 



z = (?j^(N3 + V (2,63 


®The high dimensionality of the realization is with respect to the total 
number of samples N of the process. This makes impractical a recursive 
estimation. 

®The prime • indicates the transpose. 
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The matrix x Nm] of full rank i - m is of the form 

= [0;n] (2.7) 

where n of dimension im x im accounts for a possible dynamical structure 
of the measurement system. The vector V is a (i - m) -dimensional zero mean 
independent measurement noise of covariance The output estimation error 
covariance of the vectors (2.3), indicated by P(N,i), resulting from the 
estimator % of Eq. (2.5) and the measurement system of Eqs. (2. 6-2. 7) is 

P(N,i) = (itz -K* <?i) • P(N) • (^2 - + ^ • R • (2.0J 

With the above notations, the following theorem is proven. 

Theorem 1 ; Given two processes of the form of (2.1) having covariance 
matrices (2.2) of values Pi(N) and Pa(N), respectively, and for the same 
arbitrary measurement system and estimator, as given in Eqs. (2.5) and (2.6), 
a necessary and sufficient condition for the respective estimation covariance 
matrices (2.8) to satisfy 

Pl(N,i) - P2(N,i) ^0 1 i i 5 N (2.9) 


is 


Pl(N) - P2(N) > 0 (2.10) 

Proof : Sufficient condition is immediate consequence of the quadratic 

form (2.11) 

Pl(N,i) - ?2(N,i) = A(N,i) = (<?2 " * ^(>^) * * ^0' (2.H) 

where 


A(N) = Pi(N) - P2(N) 


To prove the necessity part, expression (2.11) is rewritten as 

ea 


A(N,i) = (I! -70 • 


(?1 


• ACN) • (^21^1) 


Then it is observed that the rows of 


(^:2 


L 

'-r 


( 2 . 12 ) 


span the full (.N • m) -dimension space 


of A(N) and it is recalled that the rows of the matrix ^ are arbitrary. 


As a result of theorem 2 it follows that any model realization (1,1) 
having output covariance matrix Pm(N) greater than or equal to the process 
covariance P(N) generates, for the same arbitrary linear estimator and linear 
measurement system with additive independent noise, an upper bound of the 
process estimation error covariance (P[^(N,i) = P(N,i) s 0) in the interval 
1 £ i 5 N. 
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IVhen the measurement characteristics become available, the Kalman filter 
computed for the best^ model in a class with the previous properties is a rea- 
sonable choice for the process output estimation. Such a filter guarantees- 
the minimum (ivith reference to the chosen model) upper bound of the covariance 
of the estimation error, and it also prevents the divergence of the estimates. 
It is well knoitfn that this phenomenon may occur when, because of model approxi- 
mations, the calculated estimation covariance becomes unrealistically small 
with respect to the real one on the process (ref. 14) . 

The whole optimum filter, with reference to the entire model class, also 
could be sought. However this problem requires the computaticn of the limit of 
a sequence of models Mi in the class having covariance 

P(N) < Pj^.(N) < . . . < P^^CN) 
and this has not been considered here. 


III. GUARANTEED MODEL REALIZATION 


The computation of exact or partial realizations on the form of (1.1) 
from the output covariance ^(i,j) is essentially equivalent to the algebraic 
problem to find a decomposition for the covariance matrix 

P(N) = T(N) • A(N) • T'(N) (3.1) 


where P(N) has been defined in (2.2), 


A(N) 


rs(N-i) 


S 


■'S(0)J 


(3.2) 


~ O 

is a symmetrical block diagonal matrix with S(i) >0 of dimension mxm and 
T(N) is a nonsingular causal transformation. If T(N) is also casually 
invertible then it is called the innovation representation of y(*) and the 
independent process with covariance A the innovation process (ref. 15) . 
Recently, Akaike (refs. 16 and 17) gives of the realization problem an inter- 
pretation in terms of canonical regression analysis. His basic idea is adopted 
here and two regression schemes, the second being the extension of the first, 
are introduced to generate fixed-order realizations whose output covariance is 
greater than or equal to the process covariance P(N) (2.2). The starting 
point of the procedure is the observation that if y[^(i) is the output at 
instant i of a finite dimension model, then it can be expressed by 


^The best model measured by the trace (Pj 4 (N) - P(N)). 

^Strictly inequality because of our assumption that P(N) positive definite. 
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- 1) 


y^Ci) = KCi - 1) 


+ u (i - 1) 
res^ ^ 


y„Ci - r) 


C3.3) 


where K(i - 1) is of dimension in x r • m and r is a finite integer. The 
residual vector u^gg(i - 1), uncorrelated with the outputs prior to i - r, 
may be in general correlated with CyjJjCi - 1), , yjJj(i - r))'. With 

respect to the notation proposed in the previous section (2.2) indicate with 
P(i) and PJ^J(i) the portions of covariance matrices of process and ■"'idel in the 
interval i, i - 1, , . . ,1 and with A(i) the difference mat 

P„(i) - P(i) = A(i) (3.4) 


Introduce for the process covariance the partition 


P(iD 


■PiiCi - n ! Pi2(i - ir 

__ _ 

^ r 

_Pl2(i - 13 I P(i - 1) - 


(3.5) 


where Pn(i - 1) ='^(i,i) and Pj 2 (i - 1) = ("^^Ci. i - 1). ♦ • • »^(i»l)) are 
of dimension m x m and m x (i - l) . m. Choose the model order n = m • r^ and 
impose that the first s samples (yw(i - 1), ... , y[^|(i - s)) correlate 
with the residual. From expression (3,3) and related arguments, the model 
output covariance Pjf(i) can be decomposed in 


= 

"in, j K(i - 1) ' 
1 

• 

~S(i - 1) 1 R(i - 1)" 
1 
1 


Im j K(i - 1) ■ 
1 

_0 I I„,(i - 1)_ 


[R'(i - 1) 1 Pj^(i - 1)] 


0 1 ImCi - 1]_ 


(3.6) 


where K(i - 1) and R(i - 1) are both matrices of dimension m x (i - i) . m and 
are constrained, for all i: m(i - 1) > n to the form 


K(i - 1) = [K(i - 1) j 0] 
R(i - 1) = [R(i - 1) ! 0]. 


(3.7) 


The matrix K(i - 1) [m x n] has already been defined in Eq. (3.3), the matrix 
R(i - 1) [m X ms], ms S n accounts for the nonzero correlation terms of the 
residual vector Uj.ec(i - 1) in (3.3) with (yiw(i - 1), - • . , y]Ji(i - t))' and 
S(i - 1) [m X m] is the covariance of the residual u^g 3 (i - 1). 


^Constraining the model order to be a multiple of the number of outputs is 
imposed for notation sinqjlicity, but it is not strictly necessary. 
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The model covariance Pm(N) is therefore completely defined by a sequence 
of decompositions (3.6) or, in alternate form, with the introduction of a new 
sequence of covariance matrices V(i)[m»i x m’i], by 


V(i) 


Im+KCi-1) -R* (i-1) -S-l Ci-1) I K(i-l) 

R»(i-l)-S-^Ci-l) |ImCi-n 
iro+K(i-l) *R' (i-1) 'S-^ (i-1) ‘ K(i-l) 


S(i-l)| 0 


|V(i-l) 


and 


-1 • 


R*(i-l)*S-Hi-l) jIraCi-1) 


(3.8) 


V(i - 1) = Pj^Ci - 1) - [R» • S-1 • 


where the matrices P^Ci - 1), R(i - 1), S(i - 1) are derived from (3.6) with 
the following algorithm 


"S(i-D iR(i-l) “ 


"I„,l-K(i-1)“ 


In,l-K(i-1)' 

i 

1 

= 

t 

1 

• [Pm-R‘-S-1*R]^.j • 

1 

1 

U'(i-l)|PM(i-l)„ 


_0 ilm(i-l)_ 

LO |In,(i-l)J 


starting from 

P^CN) = V(N) = Pj^(N), and [R' • S“1 • R] = 0. 


Note also that the matrices R(i) and P[i}(i) have structure 

R(i) 

PM(i) “ PM(i) - 

for certain matrices R^(i) and Q(i). 

With the above assumptions the following theorem is proven. 

Theorem 2 ; The model representation with output covariance 

P„(N) = T(N) • A(N) * T*(N) (3.11) 

where A(N) has been defined in (3.2) with S(i) given from (3.9) and the 
transformation T(N)[N*m x N*ra] given by 


I^R(i) [m X ms] joj 
~Q(i)[m • (s-1) X m 


(s-l)]l0‘ 

I 


1 


(3.10) 


| 0 . 
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T(N) » [t(N-l)] • 


with t(i) of the form 


r>j-° _i 

Lo jt(N-2)J ' ' ' |_ 0 

KCi) • R’Cil • S“lCi)|KCij 1 

•Ci) • S-l(i) |lm • iJ 


“In, + KCi) • R’Cii 
R'Ci) • S-l(i 


guarantees an upper bound of the process output covariance matrix: 

Pm(N) - P(N) > 0 

if the matrices S(i) defined in (3.6) are chosen as 
SCO) = Ag(0) + P(l) 

. s(i) = Ag(i) + Pn(i) + [R(i) + K(i) • P^(i) - Pi2(i)l * A-ki) 

• [R(i) + K(i) • Pj^(i) - Pi2(i)]' - K(i) • R'Ci) 

- R(i) • K'(i) - K(i) • P^(i) • K'(i) V i, i = 1, . . . . N-1 

(3.12) 

for arbitrary symmetrical matrices Ag(i)[m x m] satisfying 

Ag(i) >0, i = 0, . . . , N - 2, Ag(N - 1) s 0. 

P roof : Observe first that the representation (3.11) is nothing but the 

recursive inclusions of expressions (3.8) and therefore equivalent and deduc- 
able from (3.!1>). Then decompose Pj^^CN) and P(N) according to (3.5) and (3,6) 
(i = N) . The result yields 


r[S+K-R»+R.K'+KPMK'](N-l) j (^-1)1 

L 1 Pm(n-i) J Lp12 


(N-1)|Pi2(N-1)' 


Inequality (3.13) is satisfied if 


2(N-1)|P(N-1) 


= A(N) i 0 


(3.13) 


PmCN - 1) - P(N - 1) = A(N - 1) > 0 


(3.14) 


S+K-R'+R-K»+K-P^*K'-Pn-[R+K-P|^-Pl2]*A~l*[R+K-P^-Pl2] = Ag(N-l) > 0 

(3.15) 

Proceeding backward, inequality (3.14) is successively partitioned with use of 
(3.6) and replaced by conditions of the form of (3.15) (the is replaced by 
>) . This proves the theorem. 


8 


It remains to show that (3.11) is also a causual representation of y^Ci) 
for 1 1 i 1 N and it can be reduced to a model representation of the form 
(1.1). This is easily accomplished, defining an m-valued-vector process 
(x* (i + s), x' (i + s - 1), . . . , x'(l))' with covariance matrix 


T(i + s) • A(i + s) • T* (i + s) 

= [t(i + s - D] . . . A(i + s) . . . [t(i + s - 1)]* = V(i + s). 


Then, from recursion (3.9) and (3.10), it is observed that the lowest i • m 
portion of x(*) is statistically equivalent to the corresponding yM(*) 
process from 1 to i.^® A proper restriction of x(*) with a formal shift of 
the index i of s steps on time is therefore the required realization. Note 
also that the remaining components x(j), j=i+l, . . . ,i+s are the 
predictors of the future ywC*) from observations until i according to the 
theory developed by Akaike [ref. 17). 


The model realization is finally defined by the following matrices. 
These are deducible from Eqs. (3.8) and (3.9). 


A(i) 


__ _ll 

K(i + a) 

1 

- ''1 0 « 


(3.16) 


B(i) [n X m] 


rv * K 


L UpCR' 



i+a 


C(i) 


(Oi, 



(3.17) 


(3.18) 


D(i) = Low(R» 


s-1) 


i+a 


(3.19) 


cov(w(i)) = S(i + a) 


(3.20) 


cov(xi) 



0 



loj 


1+a 


(3.21) 


^®If s = r, a noise vector has to be added to x(‘). This originates 
strictly a nonproper model realization (see eq. (3.19)), _ 

^^In the case that s • m< n - m zero matrices complete K(i) for the 
first values of i. 
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where Up(») and Low(») are respectively operators which isolate the first 
n - m and following m rows of the matrix (•), The dimension p of the 
input process w(i) results here in p = m. The index a appearing in 

(3.16)-(3,21) is a = r - 1 if s • m - n, and a = s if s • m < n - m; in 

this case, the matrices D(i) are void. 

Because of the generality of the assumptions, structures of the form of 
(3.6) can represent with a finite value of r (and s) arbitrary linear lumped 
systems. However they do not encompass the whole class of the finite realiza- 
tions, and in fact, for a given model covariance the result may not be the 

minimal realization. For stationary processes, analysis of this point has been 

at the basis of the Ho’s algorithm (ref. 1). For nonstationary processes, the 
problem received a theoretical interpretation from canonical regression 
analysis (ref. 17). Justified from the results of canonical regression analy- 
sis, an extension of our previous result is given next. 

The approach requires in this case to fix the model order n (dimension 
of the space spanned by the projection of future outputs on the past outputs) 
and the number r (m • r > n) of output samples having independent proje.-ion. 
Then at each instant i with a nonsingular transformation an n-diraensional 
basis, spanning the projection space is isolated from the whole m • r vector 
of outputs (y^Ci -!)»•••» yj^Ci - r))'. Indicate this transformation with 

rLi(i)|L2(i) j 0 1 

L(i) = i__-L (3.22) 

L 0 ImCi - 


Lj(*) [m * T X (m • r - n)], L 2 (*)[ni • r x n] . 

The model covariance assumes the following structure 


with 


Pj^(i) = L(i) • H(i) • L'(i) 


H(i) 


! HiaCi) 

"T F 

0 jHiaWjP^Ci - r)_ 


(3.23) 


(3.24) 


• r X m • r] , Hi 2 (*)[n x (i - r) • m] 

The form (structure) of H(») is explained from the fact that at every instant 
i only n independent components of the whole m • r output samples space 
correlate with the past. 

Combining the transformation (3.23) with the previous decomposition (3.6) 
results in 
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1 


~’~m 1 «i-ir 

1 

“S(i-D 1 R(i-1)“ 


-lmjK(i-l)- 

1 

.0 tL(i-l). 


1 

_R' (i-1) 1 H(i-1)„ 


_0 1 


(3.25) 

TTie recursivity of (3.25) imposes the following constraint (3.26) (this is 
derived by substituting (3.22) and (3.24) in (3.25) and developing the products) 


(k + K • Hii) • 


Lh 

CLii'kia) • Hii 


^2 • Hj 2 

ILi2 ' Hi2 
I 


= L2(i) * Hi2(i) (3.26) 


i-1 


with matrices of dimensions 

K(*)[m ^ m • r], R(*)[ni ^ m • r] 


fLiia 

LL 2 l(i 


- 1 ) I Li 2 (i - 1 ) 

a.(i - l);L2(i - 1)) =1 T 

" - n I L22(i - n 
I 


:] 


K(i - 1) = (Ki(i - l)|K2(i - D) 

Lii(*)[(ni • r - m) X (in • r - n)], Li2C*3[Cm • r - m) x n], 
L 2 i(*)[m X (m • r - n)], L 22 (*D[ni x n], 

Ki(*)[m X (m • r - n)], K 2 (*) [w ^ n]. 

The condition (3.26) is equivalent (up to a nonessential arbitrary 
transformation) 


L 2 (i) = 


K 2 (i - 1 ) 

Li 2 (i - 1 ) 


R + K • Hii 


^21 



K2(i - 1) 



* 



- 


(Lii|Li 2) • Hii 


^22 


i-1 

Li2(i - 1) 


(3.27) 


• hi2(i - 1) (3.28) 
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and 


hiaCi - l)[n x m] 


«l2(i) = [hi2(i - l)|Hi2Ci - ’)] 


(3.29) 


With the foraal substitutions of 

K • Pm • K’ I with K • H • K» I 
li li 


and 


R + K • Pm 


F 

with I I 

" L 


(R + K • Hn) « 


LI 

1 1 

1 

1- 

-- 

IK 2 • Hi2 

L2 

' J 


(3.30) 


theorem 2 is extended to realizations of the form of (3.25). From the appendix, 
where a recursion similar to (3.8) is derived for the general case, it appears 
that the matrices t(i) in (3.11) are now 


t(i) = 


Im i' K • R' • S-l 


L • R' • S' 


-1 


K • L 


m»i 


(3.31) 


From the appendix, also, the model matrices can be computed and they result in 

A(i) = In (3.32) 


B(i)[nx(m(r+i)-n)] = [(-V{ 2 *Vn ll^) • (Li 1L2)"^] I . 


i+r 


r I +K-R'*S'‘^ ! 

1 m t 

” K 

• 



[^(Ln|Li2)-R»*s-i 

t 

-Lll |Li2_ 


- vi2-v;i _ 



(i+r-i) 


C(i) = 1,22 (i + r - 1) 


(3.33) 

(3.34) 


D(i)|m X ni(r + 1) - nl 


Ch21i 1-22) 

I 


I I„, (r + 1) - n 
R' • S‘l ' 

i V' • 

I 12 


|i+r-l 

(3.35) 
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cov(w(i)) = 


[»: 


r - 1 ) I 0 

0 l^llCi + r - 1) 


covCxCD) = . V-; . 


(3.36) 


(3.37) 


IV. PROCESS APPROXIMATION 


Defined by the order n, the number m • s of the columns of the matrices 
R(*) or the number of independent output sample^ r, glasses of "guaranteed 
models" are spanned by the choice of matrices K(*)» R(*)i and Ag(*) > 0. The 
proof of theorem 2 offers a constructive sequential algorithm to generate model 
realizations. Models in these classes are general enough to encompass finite 
realizations or, in the second case, minimal realizations to every linear casual 
lumped process. Since we are interested here in the problem of the approximate 
realization, a measure of the approximation between models in the class and the 
process is defined by 


J = trace (P m(N) - F(N)) 


(4.1) 


The function J is, only for models in the class, an upper bound of the linear 
operator's norm from the Euclidean space (N • m) into itself, and is zero when 
Pm(N) = P(N) 


J > 


P„tN') 



X 

max 


(P|,,CN3 - P(N)) 


(4.2) 


where ^maxC*3 is the maximum eigenvalue of the argument matrix. With the use 
of the partitions (3.S), (3.6) and of the relation (3.14), the function J is 
given in the form 


N-l 

J = ^ ^ trace |[Ag(i) + [(R+K«P^j-Pl 2 ) (R-t-K«P^j-Pl 2 ) '3 .‘'■^sC0) | (4.3) 

i=l 


where PM(i) is generated from the recursion relationship given by (3.6) and 
u^Hi) from the following: 

A"^(i) = 




I 


[-a;i-(r+k-Pj^j-Pi2)*a 1] 


L[-Ag^*(R+K*Pj^-Pl2)‘A-i] j A-Ua-HR+K‘P^,-Pi2)’-A-1*(R+K‘Pj^,-P12)*A“^_ 

The model design is achieved by minimizing the approximation measure 

J* = min J 


li-1 


(4.4) 
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Ag(i), i = 0, N - 1 


with respect to 

R(i), i = 1. N - 1; K(i), i = I, N - 1; 

with coKstraints 

Ag(i) >0. i = 0. N - 2, Ag(N - 1) 5 0. 

A complete minimization of J is impractical, unless the value of N is 
small, A partial optimization scheme is proposed, which takes advantage of 
the special sequential structure of the performance index. For an initial 
gues£ of matrices, AsCi) ^ sequential, term-by-term, optimization with respect 
to K(i) and R(i) is performed from i=ltoi=N-l. Since the design 
matrices appear as a quadratic form on each term of the performance index, an 
analytic solution is available with use of the matrix gradient relations (ref. 
18 ). 7}^en the procedure is iterated minimizing with numerical search on 

As(i). 

Obviously, if no constraining assumptions are taken on the process, storage 
requirement and computational time limit the maximjjm size of the interval 
1, . . . , N. We note however that the proposed step-by-step optimization 
requires the inversion of matrices of fixed dimension (n x n) . The problem of 
the growing dimensionality is found on the few multiplications and sums of 
matrices required at each step to propagate P]^(i) and A"i(i). 

The possibility of extending previous results to an infinite time interval 
is investigated in the next section for a stationary process with finite 
dimension realization. 


V. APPROXIMATION WITH STATIONARY MODEL 


Introduce the partitions defined in (3.5) and (3.6) to the lowest (m + n) 
(n “ m • r) portions of process and model covariance matrices. Constrain the 
design model to be stationary. This implies that 

[S + K • R' + R • K' + K • Pm • K'IL^ = UL(PM(r)) (5.1) 

[R + K • PMlir = [UR(PMCr))jM] (5.2) 

where UL(*), UR(‘) are operators which rfspectively define the partitions of 
the upper left m x m and remaining m ' m (r - 1) upper right portions of 
the arguments and where M is an ra x m design matrix. Extend the procedure 
to i > ra + n. The conditions given in (3.12) become 

PM(r) - P(r) = A(r) > 0 (5.3) 


_^^Note that in the case of models of the foimi (3.32-3,38) we adopted R(*)» 
and K(0 must also satis^/^ the conditions (3.28). 
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UL(Pm) - Pii(r) - [(UR(P„Cr))|M) - Pi 2 (r)] * 

• [(UR(P„(r3)|M) - Pi2(r3]' = Ag(r) > 0 

C5.43 

UL(P„) - Pn(i) - [R + K • P„(i) - Pi2(i)] • 

• [R + K . Pj^Ci) - Pl2(i)]' “ ° 

(5.5) 

Note that the index i has been dropped from the matrices independent on 
It is evident from Eqs. (5. 1-5. 5) that design parameters are ACr), M, and K 
or R. Kowevex', the approximation measure is a function of the matrix A(r) 
alone, as the design matrices affect only the values of AgCO* 

J = trace |(n - r) • ULCA(r)) + A(r)| (5.6) 

the performance index is a function of ACr) alone. The other design matrices 
effect only the values of Ag(*). Therefore, the design is achieved by seeking 
for the minimumjtraeer matrix A(r), such that It is not, an empty set of values 
of M and K or R which satisfy the sequence of constraints 

Ag(i) >0, r<i<N-2, Ag(N - 1) > 0. 

Assume now to deal with stationary process having finite and stable 
realization. Expression (5,5) can be written as ULCA) - V(i) = Ag(i) where 
V(i) is given by 


V(i) = UR(A(i + 1)) • A-l(i) • UR(A(i + 1))' 

The matrix A(i) if stable ensures the existence of "guaranteed model" and is 
also a portion of the covariance of a stationary process with finite and stable 
realization. In these conditions Faurre (ref. 2) shows that sequence of 
matrices Of the form V(i) are asymptotically convergent for increasing i. 

lim V(i) = VC«) 

i-Ko 

Therefore, a sequence of approximations that are guaranteed on a finite interval 
i, . . . , N for increasing N, converges asymptotically to a guaranteed 
approximation in the infinite interval. 


VI . EXAMPLE 


To illustrate this approach, a stochastic model of the horizontal mean wind 
velocity, where altitude is taken as an independent parameter, is generated from 
the statistical analysis of historical data. The model is designed for wind 
shear prediction from onboard measurements of an aircraft on descent flight (ref. 12) 
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The wind process is defined as a two -component random vector (southerly 
and westerly) at intervals of 1 km of altitude from the ground to 10 km. Wind 
profiles are measured from twice to four times a day by launching meteorologi- 
cal balloons (rawinsonde) from the major airports. Statistical analysis of 
historical data for the areas of interest is available in the form of monthly 
or seasonal mean vectors and covariance matrices (ref. 13). A second-order 
regression scheme of the form of (3.6) with R(i) = 0 has been chosen to 
approximate a posteriori covariance of the wind components from 1 to 10 km 
conditioned to the ground measurements. Actually, the a posteriori covariance 
of the wind aloft conditioned to the ground measurement has been considered. 
The procedure proposed in section XV has been pursued constraining for simpli- 
city the sequence form 

As(i) = b(i) • Im 
and optimizing with numerical search the vector 

b = (b(0), . . . , b(N - 1))' 

The process covariance matrix and the resulting model matrices are given in 
Table 1. The normalized performance index resulted, after minimization, in 

trace (P m( 10) - P (10)) /trace (P (10)) =0.09 


VII. CONCLUSIONS 


A new approach has been proposed in this paper for the design of 
approximate realizations of stochastic processes from output covariance 
matrices in a finite interval. This approach does not require any restrictive 
assumption on the process. It appears also especially convenient in the case 
of raw experimental covariance. 

The novelty of the approach is to define, with simple sequential algorithm 
elements in classes of fixed-order dynamic models whose output covariance matrix 
bounds the process covariance in the interval of interest (the difference matrix 
is not negative definite). The design is achieved by minimizing in a chosen 
class a measure of the approximation between the model and the process. The 
specific choice of the classes, which are general enough to encompass finite 
realizations and minimal realizations of every linear process that admits 
finite dimension realization, offers two advantages: (1) the possibility to 

define a relatively simple performance measure (actually the trace of the dif- 
ference between model and process covariance matrices has been chosen) ; and 
(2) models belonging to these classes have the property that with the same out- 
put filtering or prediction structure, the estimation error covariance computed 
on the model is an upper bound of the true estimation covariance on the process. 
Apart from its intrinsic interest, when the model with such a guaranteed 
statistic is required for the simulation of the process, this property prevents 
on estimation the phenomenon of divergence often encountered when approximate 
models are used for the filter design (ref. 14) . 
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The approach, which for arbitrary process is confined on a finite-time 
interval, can furnish asymptotically convergent approximations for stationary 
finite-order realization processes on an infinite interval. 

Application is given to the modeling of mean i^rind randomness as a function 
of altitude from statistical analysis of historical wind profiles. 
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APPENDIX 


As in__the previous case of Eq, (3.8), we define a sequence of covariance 
matrices V(i), that for convenience is given here in the following transformed 
form: V(i) = U{i) • V(i) • U'(i3, with 


V(i) 



Vii{ 

0 

._L 0 

a 

0 

• VI 2 • Vil . Vj,| Hi2 


0 1 
L 1 

HI 2 

1 


^(m.r-n) j °"| 

0(1) = “--zrr- 

1^2 • 


Vii(‘)[Cin • r - n) X (m • r - n)], Vi 2 (*)[(>n • r - n) x n], V 22 [n x n], 

A recursion equivalent and deducible from (3.25) is 


L(i) ♦ U(i) • V(i) . U*(i) . L'(i) = 


flm + K • R» • S-i|k • L-n 

I 0 I L • U • V • U’ • L’ J 

■[- 


and 


[U • V • U‘] 


Im + K • R» • S'^ 1 K • 

_l 


=HCi-l) - [S' • 


1-1 


where the matrices H(*), R(*), S(») are derived from (3.25) with the following 
algorithm 


r S I R • L' 1 flm I -K • L-n 

— 7T-7-“ = -f-— . [L* 

I L • R* jL • H • L* . , 0 I If. ... . , 

L ' -J 1-1 •- I m(i-l)_i 1-1 


(H - R’ • • R) • L'3 


E 


Im ! -K • I-*' 


° J ^m(i-l) J |i-l 
18 


starting from 


L(N) • H(N) • L'(N) = Pj^(N) and [R» • • R] = 0. 
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TABLE l.-WIND COVARIANCE - MODEL MATRICES 


A Poitonoii Covariance Matrix of the Two Wind Componenlt (Weitorlv and Soulherly) from 10 to 1 km of Allilude, 

Condillonad to the Ground Meaiurenwmi 


km 


CAPE KENNEDV, FLORIDA 
JANUARY 


PERIOD OF DATA 

Jan. 1, 1950 to Nov, 17, 1956 
Nov. 18, 1956 to Dec. 31. 1963 


10 


0.237E 

02 



















0.232E 

Q'i 

0.207E 

02 

















O.tOSE 

02 

O.308E 

01 

0.318E 

02 















0.229E 

01 

0.206E 

02 

0.572E 

01 

0.3S8E 

02 













0.172E 

02 

0.401E 

01 

0.207E 

02 

0.893E 

01 

0.40DE 

02 











0.6Q5E 

00 

0.209E 

02 

0.3B2E 

01 

0.37BE 

02 

0.7B6E 

01 

0.6246 

02 









0.15BE 

02 

0.431E 

01 

0.207E 

02 

0.104E 

02 

0.370E 

02 

0.104E 

02 

0.4S3E 

02 







-0.920E 

00 

0.212E 

02 

0.268E 

01 

0.3B1E 

02 

0.804E 

01 

0.6U9E 

02 

0.102E 

02 

0,6196 

02 





0.160E 

02 

0.377E 

01 

0.2S0E 

02 

0.1 10E 

02 

0.378E 

02 

0.1276 

02 

0.4G7E 

02 

0.142E 

02 

0.604E 

02 



-O.210E 

01 

0.227E 

02 

0.170E 

01 

0.377E 

02 

0.608E 

01 

0.62GE 

02 

0.10BE 

02 

0.628E 

02 

0.156E 

02 

0.759E 

02 

0.149E 

02 

0.407E 

01 

0.2SBE 

02 

0.142E 

02 

0.382E 

02 

0.1B3E 

02 

0.4676 

02 

0.174E 

02 

0.6046 

02 

0.201E 

02 

0.738E 

02 



















-0.203E 

01 

0.240E 

02 

0.808E 

00 

0.403E 

02 

0.481E 

01 

0.6S0E 

02 

0.S68E 

01 

0.643E 

02 

0.1446 

02 

0.7736 

02 

0.107E 

02 

0.93BE 

02 

















0.1B3E 

02 

O.OGOE 

01 

0.297E 

02 

0.172E 

02 

0.409E 

02 

0.107G 

02 

04976 

02 

0.2306 

02 

g.63GE 

02 

0.265E 

02 

0.772E 

02 

0.281E 

02 

0.9SOE 

02 















-0.232E 

01 

0.23BE 

02 

0.230E 

01 

0.389E 

02 

0.7276 

01 

O.S806 

02 

0.108E 

02 

0.64BE 

02 

0.1826 

02 

0.7846 

02 

0.223E 

02 

0.944E 

02 

0.31BE 

02 

0.111E 

03 













. 0.14BF 

02 

0.013E 

01 

0.297E 

02 

0.201E 

02 

0.4196 

02 

0.246E 

02 

0.492E 

02 

0.281E 

02 

Q.623E 

02 

0.31 2E 

02 

0.76SE 

02 

0.324E 

02 

0.96GE 

02 

0.379E 

02 

0.11BE 

03 











-0.478E 

01 

0.237E 

02 

0.143E 

01 

0.410E 

02 

0.737E 

01 

0.B82E 

02 

0.111E 

02 

0.670E 

02 

0.186E 

02 

0.8176 

02 

0.237E 

02 

0.074E 

02 

0.329E 

02 

0.1 12E 

03 

0.40BE 

02 

0.129E 

03 









0.150E 

02 

0.788E 

01 

0.319E 

02 

0.230E 

02 

0.44BE 

02 

0.2966 

02 

0.5116 

02 

0.3206 

02 

0.642E 

02 

0,371 E 

02 

0.772E 

02 

0.3B4E 

02 

0.996E 

02 

0.460E 

02 

0.1 20E 

03 

0.E08E 

02 

0.147E 

03 







-0.588E 

01 

0.232E 

02 

0.013E 

00 

0.398E 

02 

0.717E 

01 

0.579E 

02 

0.966E 

01 

0.679E 

02 

0.177E 

02 

0.842E 

02 

0.231E 

02 

0.99DE 

02 

0.32SE 

02 

0.1 12E 

03 

0.4156 

02 

0.1306 

03 

0.505 E 

02 

0.1506 

03 





0.141E 

02 

0.732E 

01 

0.298E 

02 

0.2B7E 

02 

0.442E 

02 

0.31BE 

02 

0.5216 

02 

0.391 E 

02 

0.666E 

02 

0.449E 

02 

0.777E 

02 

0.490E 

02 

0.991 E 

02 

0.S47E 

02 

0.116E 

03 

0.6CBE 

02 

0.1BOE 

03 

0.619E 

02 

0.182E 

03 



-0.71SE 

01 

0.204E 

02 -0.1 39E 

01 

0.3B1E 

02 

0.9B2E 

01 

0.530E 

02 

0.7766 

01 

0.6356 

02 

0.1716 

02 

0.7956 

02 

0.227E 

02 

0.948E 

02 

0.327E 

02 

0.103E 

03 

0.4266 

02 

0.12SE 

03 

0.504E 

02 

0.14BE 

03 

0.632E 

02 

0.163E 

03 


Modal Matricoi from 10 to 1 km of Allituda 
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CO) - 1 
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«>0B tjmUTY 


10 

0 

8 

7 

6 

5 

4 

3 

2 


AO) 


cov IwO)) 


0.7860 00 

0.45204)2 

0.3440 

02 

0.1160 

01 

0.5250411 

0.820D 00 

0.1160 

01 

0.3460 

02 

0.7810 00 

0.17104)1 

0.2361? 

'.)2 

-0.108D 

01 

0.62704)1 

0.8210 00 

-0.1080 

01 

0.2370 

02 

0.8230 00 

-0.9730-02 

0.2140 

02 

0.4830 

00 

0.3100411 

0.8410 00 

0.4830 

00 

0.1840 

02 

0.7910 00 

-0.24004)1 

0.1510 

02 

-0.1280 

00 

0.787002 

0.8570 00 

-0.1280 

00 

0.1880 

02 

0.8280 00 

-0.12904)1 

0.1480 

02 

0.2220 

00 

0.56004)1 

0.8180 00 

0.2220 

00 

0.1410 

02 

0.7980 00 

-0.48304)1 

0.1110 

02 

0.1020 

01 

0.42404)1 

0.7990 00 

0.102D 

01 

0.122b 

02 

0.8640 00 

-0.2930-01 

0.1360 

02 

0.7390 

00 

0.4050411 

0.8200 00 

0.739D 

00 

0.1610 

02 

0.7360 00 

-0.6830-01 

0.1300 

02 

0.2190 

00 

0.72204)1 

0.6770 00 

0.219D 

00 

0.1220 

02 

0.6050 00 

-0.1140 OO 

0.1940 

02 

0.342D 

01 

-0.1240 00 

0.5820 00 

0.3420 

01 

0.1860 

02 


«>»(X(|0)) 


0.184E 03 0.632E 02 
0.632E 02 0.165E 03 
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